A hierarchy of self-consistent stochastic boundary 
conditions for Ising lattice simulations 



Yidan Wang 3 , You Quan Chong a b , Siew Ann Cheong a * 

"Division of Physics and Applied Physics, School of Physical and Mathematical Sciences, 
Nanyang Technological University, 21 Nanyang Link, Singapore 637371, Republic of Singapore 
b Physics and Astronomy Department, Stony Brook University, Stony Brook, NY 1 1794-3800, 

United States of America 



Abstract 

We describe a hierarchy of stochastic boundary conditions (SBCs) that can be 
used to systematically eliminate finite size effects in Monte Carlo simulations of 
Ising lattices. For an Ising model on a 100 x 100 square lattice, we measured the 
specific heat, the magnetic susceptibility, and the spin-spin correlation using SBCs 
of the two lowest orders, to show that they compare favourably against periodic 
boundary conditions (PBC) simulations and analytical results. To demonstrate 
how versatile the SBCs are, we then simulated an Ising lattice with a magne- 
tized boundary, and another with an open boundary, measuring the magnetization, 
magnetic susceptibility, and longitudinal and transverse spin-spin correlations as 
a function of distance from the boundary. 
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1. Introduction 

We have learnt a great deal about the statistical physics of critical phenomena 
from Monte Carlo simulations of the Ising model. Starting with early works 
(see Ref. i4\ for a review), Monte Carlo simulations have helped us verified exact 
solutions for the two-dimensional Ising model o, |6(] (see Ref. []7|] for a review) 
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as well as renormalization-group calculations for higher-dimensional Ising mod- 
els AH] (see Ref. [10] for a review). Associated with the growing interest in 
Monte Carlo methods within the statistical and computational physics communi- 
ties, there were also notable developments of important algorithms to tackle the 
problem of critical slowing down, where the dynamical time scale diverges as we 
approach the critical point [| 1 lM 1 3fl . As simple Monte Carlo algorithms become 
highly inefficient, many algorithmic acceleration methods based on resampling or 
cluster moves have been proposed |[l4l - l27ll . A good review of the basic principles 
behind acceleration algorithms can be found in Ref. |28ll. Recently, there have 
also been developments in parallelization acceleration Il29ll and hardware acceler- 
ation [3(^-32|]. These represent improvements beyond the previous state of the art 
in Monte Carlo simulation of Ising models (1331 - 43511 . 

The origin of this diverging dynamical time scale scale at the critical tem- 
perature is the diverging correlation length scale [36, [13]. Existing acceleration 
algorithms addresses only the problem of the diverging time scale. To address the 
problem of the diverging length scale, we either simulate a system whose size is 
larger than the diverging length scale, or perform finite size scaling analysis to 
eliminate finite size effects I138l - l44l] . In nearly all cases, we follow suggestions in 
texts on Monte Carlo methods to impose periodic boundary conditions 

(PBCs). In other studies, artificial correlations introduced by PBCs have been 
observed to strongly influence the results of computational studies I48"l-l50ll. 

In this paper, we propose a hierarchy of self-consistent stochastic boundary 
conditions (SBCs) for the Monte Carlo simulation of Ising lattices. The idea of 
designing boundary conditions that allows us to minimize, or completely elim- 
inate finite size effects is not new. For exact diagonalization studies or quan- 
tum Monte Carlo simulations, twist boundary conditions have been introduced 
11591 - 16 ill . In the Monte Carlo literature, Binder and Muller-Krumbhaar first used 
self-consistent boundary conditions in their study of a simple classical Heisen- 
berg ferromagnet HH). In this boundary condition, a uniform magnetic field is 
applied to the boundary, and iteratively equilibriated with the average magneti- 
zation of the interior spins. Hasenbusch then introduced fluctuating boundary 
conditions, for his Monte Carlo study of the three-dimensional Ising model 115 ill , 
by averaging over periodic and antiperiodic boundary conditions. These fluctuat- 
ing boundary conditions were generalised [1521, 15311. by introducing phase shifts in 
boundary spins of quantum lattices. These boundary phase shifts are adjusted it- 
eratively, until measurements at the boundaries are consistent with measurements 
within the bulk. In addition, the name self-consistent boundary conditions has 
also been used by Olsson [1551. 15611 to describe his method of interpolating between 
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periodic and fluctuating boundary conditions to achieve self-consistency between 
the boundary and the bulk in simulations of the XY model. Our contributions in 
this paper is the development of boundary conditions that are stochastic and self- 
consistent, and with no freedom to perform optimization. Finite size effects can 
be systematically eliminated by going to higher and higher order in the hierarchy 
of stochastic boundary conditions. More importantly, our SBCs can be used to 
simulate asymmetric lattices with one or more special boundaries. 

Our paper is organized as follows. In Section [21 we describe how our hier- 
archy of SBCs can be generated by sampling the spin flip statistics of larger and 
larger clusters, and how we re-sample flips of the boundary pseudospins from 
these distributions. We then describe how we refresh the spin flip statistics so 
that the simulation becomes self-consistent eventually, which we check through 
simple measurements. In Section [31 we compare the performances of the zeroth- 
order (SBCO) and first-order (SBC1) SBCs against that of the PBC, by measuring 
the specific heat, magnetic susceptibility, and spin-spin correlation function of the 
Ising lattice. In Sectional we simulate two asymmetric lattices using SBCO and 
SBC1. In the magnetized boundary simulation, one boundary of the Ising lattice 
is coupled to pseudospins which are always up, to simulate the effect of a short- 
range DC magnetic field at the boundary. In the open boundary simulation, one 
boundary of the Ising lattice is left free, i.e. not coupled to any pseudospin, to sim- 
ulate an actual surface. Along the other three boundaries, we impose SBCs with 
spin flip statistics that can self-consistently vary with distance from the special 
boundary in both cases. We then conclude in Section [5J 

2. Methods 

2.1. Overview of Boundary Conditions 

In a computer simulation, the number of spins N is necessarily finite. Unlike 
the infinite Ising lattice considered in the thermodynamic limit, a finite lattice of 
spins presents boundaries, which must be treated with care. If we choose to sim- 
ulate an Ising lattice using open boundary conditions, then as shown in Figure 
Etta), different spins will have different coordination numbers. The translational 
symmetry of the infinite Ising lattice will thus not be preserved. In practice, PBCs 
are imposed, as suggested by textbooks on Monte Carlo methods BH. With 
PBCs, all spins have the same coordination number, as shown in Figure [TJb), and 
are thus translationally equivalent. However, spins on the boundaries become ar- 
tificially correlated with spins on the opposite boundaries. When the correlation 
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length is much shorter than the size of the simulation system, these artificial cor- 
relations do not affect measurements performed in the Monte Carlo simulation. 
However, they will pose a problem when the correlation length becomes com- 
parable to the size of the simulation system, at temperatures close to the critical 
temperature T c . Typically, finite size scaling analyses are performed to extract the 
infinite-lattice limits of measured quantities OSl - Ml . 




(b) (c) 



Figure 1 : Simulating a 4x4 Ising lattice using (a) open boundary conditions; (b) periodic boundary 
conditions; and (c) stochastic boundary conditions. With open boundary conditions (a), the coor- 
dination numbers, three and two respectively, of the boundary spins and corner spins are smaller 
than the coordination number, four, of the interior spins. With periodic boundary conditions (b), 
all spins are coupled to four other spins, but spins on the boundaries will be artificially correlated 
with spins on the opposite boundaries. These artificial correlations are removed in the stochastic 
boundary conditions (c), where boundary spins are coupled to independent pseudospins. 

If instead of PBCs, we couple each spin on the boundaries to an independent 
pseudospin, as shown in Figure (He), it is also possible to ensure that all system 
spins are coupled to the same number of nearest neighbors. Pseudospins do not 
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contribute towards measurements, which are done only over the system spins, but 
must also be updated from time to time, to simulate the coupling of the system 
of N spins to a fluctuating environment. Our goal is to choose an appropriate 
pseudospin dynamics to mimic the fluctuating environment seen by a subsystem 
of N spins within an infinite Ising lattice at temperature T. 

2.2. Hierarchy of Stochastic Boundary Conditions 

Clearly, if we could simulate an infinite Ising lattice, we would simply extract 
an ensemble of histories of the environment spins coupled to the Af-spin system, 
as shown in Figured and use these as our SBCs. Such a set of SBC would in fact 
be exact, i.e. the values of all quantities that can be measured entirely within the 
A^-spin system would be identical to their infinite-lattice values. This is because 
all possible correlations between the A^-spin system and its infinite environment 
have been incorporated into the dynamics of the SBCs. 




Figure 2: From Monte Carlo simulations of the infinite Ising lattice, the dynamical histories of the 
16 environment spins (colored light gray) can be used as stochastic boundary conditions for the 
4x4 system they are coupled to. 

In practice, the SBCs could not be devised in this manner. We should also 
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not simulate a larger Ising lattice (with PBCs), within which the Af-spin system 
is embedded, extract the dynamical histories of the relevant environment spins, 
and use these as our SBCs. These SBCs capture exactly the correlations between 
the iV-spin system and its finite environment, and therefore will not do better than 
direct simulations of the larger finite lattice in approximating the infinite-lattice 
limits of various observables. 

However, by observing translational and rotational symmetries of the infinite 
Ising lattice, we can develop a hierarchy of SBCs that approximates the infinite- 
lattice embedding better and better. First, let us observe that the infinite Ising 
lattice is translationally invariant, and thus the dynamics of the environment spins 
must be statistically identical to the dynamics of the system spins. At the zeroth 
order, an environment spin must flip no faster, or no slower than a system spin. If 
we call the time interval between two consecutive flips (T-to-J, or J,-to-T) of a given 
spin the flip time of the said spin, the flip time distributions of an environment 
spin must be identical to that of a system spin. Therefore, to mimic the fluctuating 
environment of a iV-spin system we design an SBC whereby the pseudospin flip 
statistics are identical to the system spin flip statistics. We call this a zeroth- 
order stochastic boundary conditions (SBCO), because no spin-spin correlations 
are taken into consideration. 

In general, for the ferromagnetic Ising model, a I spin surrounded by T spins 
is much more likely to flip compared to an f spin surrounded by | spins. In 
contrast to the plain spin flip statistics, the spin flip statistics conditioned by the 
spin configuration of its nearest neighbors contain information on the spin- spin 
correlations. Therefore, we would like to flip a pseudospin more frequently if it 
is misaligned with the system spin it is coupled to, and less frequently otherwise. 
Again, because of translational symmetry in the Ising model, the pseudo-system 
spin pair should have the same dynamics as any nearest-neighbor pairs within the 
system. If we design the SBCs such that the conditional pseudospin flip statistics 
are identical to the conditional spin flip statistics, we expect to capture part of 
the spin-spin correlations in the infinite Ising lattice. We call this a. first-order 
stochastic boundary conditions (SBC1), because correlations between nearest- 
neighbor spins are partly accounted for. 

Going further, we note that next-nearest neighbor spins make more and more 
significant contributions to the correlations between spins as the correlation length 
increases. To capture this correlation, we need to go to higher-order conditional 
spin flip statistics, whereby the spin configuration of next-nearest neighbors, in 
addition to nearest neighbors, are taken into consideration. Each pseudospin has 
one nearest-neighbor system spin, and two next-nearest-neighbor system spins. In 
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principle, the pseudospin flip statistics must be conditioned on the configuration 
of these three spins, for which there are eight. Once again, these conditional pseu- 
dospin flip statistics should be identical to those measured within the system. We 
call this the second-order stochastic boundary conditions (SBC2). We expect to 
capture more of the spin-spin correlations within the infinite Ising lattice, because 
correlations between next-nearest-neighbor spins are also partly accounted for. 




Qth j st 2 nd 3 rd 4 th 



Figure 3: Hierarchy of spin clusters whose conditional spin flip statistics can be applied to stochas- 
tic boundary conditions. 

In this way, a hierarchy of SBCs that will better and better approximate a 
finite system of spins embedded within an infinite environment can be built up, as 
shown in Figure [3] At some point, our measurements will so closely approximate 
the infinite-system results that we do not refine the SBCs any further. At high 
temperatures, where the spin-spin correlation length is short, we will show in 
Section [3] that we do sufficiently well with zeroth and first order SBCs. Close 
to T c , we will need higher-order SBCs. In addition to finite size scaling, where 
simulation systems of different sizes are employed, we can also do finite order 
scaling, where we simulate finite Ising lattices using SBCs of different orders, and 
then extrapolate to infinite order. 

2.3. Overview of Algorithm 

Unlike for PBC, we can only start a SBC simulation after first knowing the 
conditional spin flip statistics. These can be obtained from theory, but we pre- 
fer not to. Instead, we will measure the spin flip statistics in a calibration stage, 
where the finite Ising lattice is simulated using PBC. After the system equilibri- 
ates, i.e. total energy and magnetization become stable, we start collecting flip 
time statistics. In this calibration stage, we will need a long time to accumulate 
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enough statistics if we monitor only a single observation cluster, which typically 
flips after a large number of Monte Carlo moves. Since the infinite Ising lattice is 
translationally invariant, spin flip statistics from observation clusters at different 
sites should be identical. Furthermore, since the Ising lattice is rotationally invari- 
ant, spin flip statistics from observation clusters at different orientations should 
also be identical. We can therefore take advantage of these symmetries to ac- 
cumulate spin flip statistics faster from multiple (and overlapping) observation 
clusters. 

Once the spin flip statistics have been measured, we will turn off PBC, and 
impose SBC. We call this the implementation stage. For a L x L square Ising 
lattice coupled to 4L pseudospins, the Hamiltonian is given by 



Here, = ((i u i 2 ), (71,72)) denotes nearest-neighbor system sites, whereas 

(k, k') = ((ki,k 2 ),(k' l ,k' 2 )) denotes nearest-neighbor system-pseudospin sites. Sj = 
+ 1 represents an | spin while Si = -1 represents a I spin at site i (pseudospin site 
is indicated with a prime). The Metropolis algorithm is used to flip the system 
spins, while the pseudospins are updated using the SBC algorithms described be- 
low. 

2.3.1. Zeroth-Order SBC 

For SBCO, we need only measure the unconditional spin flip statistics. Moni- 
toring the spin at i for a long time, we will find it flipping at times (t\,t 2 , . . . , t M )- 
This time series is autocorrelated, more strongly so when the system is closer to 
T c . In SBCO, we will assume the spin flips to be uncorrelated, and determine the 
flip times {At\, . . . , At M -\), where At m = t m+ \ - t m . We then store the f-to-J, flip 
times in queue A and the i-to-t flip times in queue B. 

After enough flip time data is collected in the calibration stage, we can impose 
SBCO and start the implementation stage, by randomly initializing each pseu- 
dospin to be t or J,. If a pseudospin k' is T, we then randomly draw one waiting 
time At from queue A. Else, we draw At from queue B. We will then wait At time 
steps, before flipping the pseudospin k', and draw a new waiting time At from the 
appropriate queue. 

To achieve self-consistency, both queues are implemented in C as arrays with 
fixed length M. This data structure mimics a first-in-first-out queue for data collec- 
tion, at the same time offering us the advantage of random access for resampling 
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waiting times. After the queues are filled up from the front, we continue to ac- 
quire spin flip statistics from within the system. The newest flip time data that we 
obtain each time will overwrite the oldest flip time data with the aid of a moving 
index, as shown in Figure |4] In this way, only the M most recent flip times are 
stored in each queue, to ensure that the simulation 'forgets' that it started with 
PBC. Measurements are taken only after the simulation becomes self-consistent. 




Figure 4: Overwriting the oldest flip time with the newest flip time with the aid of a moving index 
in a fixed-length array. 



2.3.2. First-Order SBC 

For SBC1, instead of single spins, we consider spin pairs and distinguish be- 
tween four different cases. For a given spin pair 5i5 2 , we call the first spin the 
target spin, and the second spin its neighbor spin. We measure the spin flip statis- 
tics of the target spin, conditioned on the state of the neighbor spin. To store the 
flip time data of the target spin, we now need four different queues, A T (Tt-to- 

XT)» #t (XT-to-TT). A± (TX-to-XX), and B i (XX-to-TX)- Because there are now four 
instead of two queues to fill, we must collect more spin flip statistics from the sim- 
ulation. Therefore, whenever a system spin is flipped, we use all four spin pairs it 
is part of to update the queues. 

To understand how the conditional spin flip statistics is collected for SBC1, 
consider the two examples shown in Figure [H In Figure [5|a), the left target spin 
flips from f to J, at time t\, and then back to t at time t 2 . Throughout this time 
interval, the neighbor spin remains f ■ Therefore, we store the flip time At = t 2 - 1\ 
in the queue Bp In FigureOb), the neighbor spin flips from T to J, at time t\ < t 3 < 
t%. Therefore, the |-to-t spin flip of the target spin at t% is conditioned tf = £3 - t\ 
of the time by a t neighbor spin, and = t% - h of the time by a J, neighbor spin. 
We add the flip time At = t% - t\ to two queues, B^ with weight Wf = m At, and 
Bi with weight wj, = tJAt. We can also understand the data collection for Figure 
Oa) in terms of these weights: since t± = 0, Wf = 1 and = 0. 

As with SBCO, pseudospins will be randomly initialized at the start of the 
implementation stage. There are then two ways to draw a waiting time for a pse- 
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Figure 5: Measuring the flip time of the left target spin when (a) the right neighbor spin stays 
spin-t between t\ and ti\ and (b) when the right neighbor spin flips at some t\ < £3 < ti. 



duospin. In the first approach, we first draw a continuous random number W from 
the uniform distribution U(0, W M ), where W M = Z^i Wt lS the total weight stored 
in the flip time array. We then run through the weight array until the cumulative 
weight Yl t =\ w t ^ W. The f*th entry in the flip time array will be our new waiting 
time. Alternatively, we can start by selecting a random entry in the flip time array, 
whose weight is w. We then draw a continuous random number r from U{0, 1). 
If r < w, we accept this entry as our new waiting time. Else, we will select a 
new random entry and a new U(0, 1) random number r, until the random entry is 
accepted. In both approaches, the probability of a waiting time being selected is 
proportional to its weight in the flip time array. The first approach uses the random 
number generator efficiently, because every trial is accepted, but is slow because 
of the need to sum the weights. The second approach wastes some random num- 
bers, because not every trial is accepted, but is fast because no cumulative sums 
are evaluated. In this study, we adopted the second approach for generating wait- 
ing times. 

More importantly, while we are waiting for a f pseudospin to flip, its neighbor 
system spin may flip from f to I, or vice versa. If the neighbor system spin had 
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remained f (J,) throughout, we should sample the pseudospin waiting time A? TT 
(A? T j_) from A T (A|). When the neighbor system spin flips back and forth between 
t and i, the proper waiting time A/f should be a linear combination of A% and 
Ar TJ .. We also expect the contribution from At-^ (At^) to be proportional to the 
total time the neighbor system spin is | (J,) during the waiting time. However, 
since we do not know when the neighbor system spin will flip, we cannot calculate 
the waiting time At^ right after the pseudospin flipped. 

To solve this problem, we introduce a waiting fraction f w , which is the fraction 
of total waiting time completed by a pseudospin. Right after a pseudospin flipped, 
we draw A% and At n from A T and Aj,, and set f w = 0. For subsequent time 
steps, if the neighbor system spin is J, f w increases by 1/A? TT , whereas if the 
neighbor system spin is |, f w will increase by 1 /At^. In this way, f w will increase 
at an average rate determined by the proportions of times the neighbor system spin 
spends in the f and J, states. When f w reaches or exceeds 1, we flip the pseudospin. 

2.4. Self-Consistency 

At the start of the SBC simulations, the boundary conditions are inconsistent, 
because we are flipping pseudospins using spin flip statistics obtained from PBC. 
These would be contaminated with artificial correlations. Therefore, while the 
SBC simulations are running, we keep measuring the spin flip statistics, and use 
the updated spin flip statistics for the SBC simulations. Because spins on opposite 
boundaries are no longer coupled to each other, artificial PBC correlations will 
die out, and we will be left with whatever approximation of the infinite-lattice 
correlations the SBC admits. Eventually, the spin flip statistics will no longer 
change with further updating, and we say that the SBC simulation has become 
self-consistent. Measurements can then begin. 

To check that artificial PBC correlations in the flip time statistics collected 
from the calibration stage indeed die out in the implementation stage, and how 
long it takes for the system to become self-consistent, we examine the distributions 
K{t) of flip times in the data queues. We measure changes in these distributions 
in two ways, by (1) computing the Jensen-Shannon divergence (JSD) between 
the current distribution and the distribution at an earlier time, as well as (2) by 
measuring the moments of the old and new distributions. 

2.5. Measurements 

Once the SBC simulation is self-consistent, we start the measurements of var- 
ious physical quantities. We do this in the high temperature regime, where low 
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order SBCs are expected to be able to capture most of the weak correlations be- 
tween spins, for three scenarios: 

1 . an L x L Ising lattice embedded within an infinite system. Such simulations 
are labeled SBCO and SBC1 depending on whether the zeroth-order or first- 
order algorithm has been used; 

2. an L x L Ising lattice embedded within a semi-infinite system, such that the 
pseudospins along one boundary are perfectly magnetized. The other three 
boundaries are coupled to a fluctuating environment. Such simulations are 
labeled SBCO + M and SBC1 + M depending on which SBC algorithm has 
been used; and 

3. an L x L Ising lattice embedded within a semi-infinite system, such that one 
boundary is completely open (i.e. not coupled to pseudospins). The other 
three boundaries are coupled to a fluctuating environment. Such simulations 
are labeled SBCO + O and SBC1 + O depending on which SBC algorithm 
has been used. 

The purpose of the first scenario is to compare the performances of the SBC 
simulations against the PBC simulation, as well as against analytical results. We 
simulate the second and third scenarios because these cannot be easily done using 
PBC, to demonstrate the versatility of SBCs. 

2.5.1. Correlation Time 

The first quantity we measure in each simulation is the correlation time r. In 
each Monte Carlo time step, at most one spin is flipped and physical properties 
of the system remain largely the same. Therefore the data collected for consecu- 
tive time steps will be highly correlated. To obtain statistically independent data 
points, we need to let the system evolve for a time on the order of t. 

To measure r, we compute the autocorrelation function 

j 'max - ' 

X {t) = — Yj m(f)m(t' + t) - <m> 2 (2) 

of the average magnetization m, and fit it to a decaying exponential of the form 

X (t)=x(0)^p(~)- (3) 

For the rest of the physical quantities, we then perform statistically independent 
measurements every 2t 1450. 
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2.5.2. Other Measured Quantities 

For the first scenario, where our goal is to compare SBC simulations against 
PBC simulations, we measure the specific heat 



c(T)=/3 2 N((e 2 )-(e) 2 ), 



(4) 



and the magnetic susceptibility 



X(T)=/3N((m 2 )-(m) 2 ), 



(5) 



which are defined in terms of the variances of the average energy per spin 



Here, /3 = 1 /k B T, with the Boltzmann constant k B = 1 for our simulations, T 
is the temperature, and N is the total number of system spins. The notation (i, j) 
indicates that i and j are nearest neighbors while (k, k') indicates the pseudospin at 
k' is a nearest neighbor of the spin at k. Finally, (•) indicates an ensemble average 
over statistically independent samples from multiple simulations, 
We also measure the spin-spin correlation 



between two spins at i and j. Since the Ising lattice is translationally and rotation- 
ally invariant, we can average G(i, j) over all pairs (i, j) with the same separation 
r = |i-j| to get 




(6) 



and the average magnetization per spin 




(7) 



G(i,j) = (5 i 5 j )-<5 i )(5 j ). 



(8) 




(9) 



where n(r) is the number of pairs with separation r. 
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2.5.3. Measurements for Magnetized and Open Boundaries 

For the magnetized and open boundaries, we would like to see how a single 
special boundary affect physical quantities like the magnetization per spin and 
magnetic susceptibility at different distances from the boundaries. To do this, we 
measure the magnetization per spin 



i rows away from the special boundary. The average magnetization and magnetic 
susceptibility i rows away from the special boundary are then (ra(z')) and 



Here L is the number of spins in one row. We do this one row at a time because the 
presence of a special boundary breaks the translational invariance in the direction 
perpendicular to the boundary, but translational invariance is maintained in the 
other direction. 

We also expect the spin-spin correlation to vary from row to row because of 
the special boundary. More importantly, rotational invariance is also broken by the 
presence of the special boundary, so we have two distinct limits for the spin-spin 
correlation function. We call correlations between spins in the same row (same 
distance from the special boundary) the transverse spin-spin correlation G ± (r) 
and correlations between spins in the same column (different distances from the 
special boundary) the longitudinal spin-spin correlation Gy(r). 

3. Performance of SBC 

3.1. Self-Consistency 

For SBCO simulations, each datum has the same weight. If queue A has M 
entries, and the flip time t appears n{t) times, then K{t) = n(t)/M. For SBC1 
simulations, each datum has a different weight. The total weight of the data in 
queue A T is = ^tOX where w T (i) is the weight of the z'th flip time stored in 
the queue A T . If the flip time t occurs in entries t\,t 2 ,... t n{t) , the probability K^(t) 
for finding the target spin flip from j to I in t time steps given its neighbor spin is 
t will be given by 




(10) 




(11) 




(12) 
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To test for self-consistency, we measure the distributions K n (t) at a set of dis- 
crete times {t„}. The time interval At = t„ + i - t„ is chosen to be large enough for 
a significant fraction of the queues to be refreshed, but small enough that we can 
still monitor the long-time evolution of K n (f). 

3.1.1. Jensen-Shannon Divergence 

For two flip time distributions K\(f) and K 2 (t), where t = 1,2, ... , we can 
define their Jensen-Shannon divergence (JSD) to be 

OO ^ CO - CO 

JSD(^,^ 2 ) = -J^K(t)\nK(t) + -J^K l (t)\nK l (t) + -J^K 2 (t)\nK 2 (t). (13) 

r=i t=i t=\ 

Here, 

K(t) = ^[K 1 (t) + K 2 (t)] (14) 

is the average flip time distribution obtained by combining the flip time statistics 
for Ky(f) and K 2 (f). The Jensen-Shannon divergence [62], which is a symmetrized 
version of the Kullback-Liebler divergence rf63l l64ll. measures the statistical dis- 
tance between two probability distributions. The larger the JSD, the more different 
the two distributions are. 

In our study, we calculate JSD of successive flip time distributions K n (t) and 
K n+ i(t) to observe possible changes in the distributions. We expect 

JSD„ = JSD^,^) (15) 

to decrease initially, and fluctuate about a constant value after the flip time dis- 
tribution has converged. However, when we plot the JSD value of the T-to-J, flip 
time distributions for a SBCO simulation of a 100 x 100 Ising lattice at T = 3.0 as a 
function of time in Figure [6l we see that the JSD value fluctuates weakly about an 
average value of JSD = 0.0404. To better appreciate this average JSD, let us note 
that the numerical value of the JSD represents the difference between two distri- 
butions in number of bits. There is thus about 0.04 bits of difference, or 4 x 10~ 8 
bits of difference per data points, between successive flip time distributions, even 
though a significant fraction of the queue has been refreshed in the intervening 
time. In particular, this suggests that the flip time distribution estimated from the 
PBC calibration stage is already very close to the self-consistent flip time distri- 
bution. To confirm this, we also calculated the moments of the distributions. 
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Figure 6: The JSD value for the T-to-J, flip time distributions versus the number of sweeps of a 
SBC0 simulation of a 100 X 100 Ising lattice at T — 3.0. The length of the queue used to collect 
the T-to-J, flip time data is M — 1, 000, 000, and the time interval between successive distributions 
is At = 200 sweeps. Each sweep consists of 10,000 Monte Carlo steps. 

3.1.2. Moments of flip time distributions 

The kih. moment fi k of a flip time distribution K(t) is 

lu k = J]t k K(t). (16) 

t 

Any change in K(t) would therefore be reflected as changes in its moments. In- 
specting the first four moments for the T-to-J, flip time distribution at T = 3.0 for 
a 100 x 100 lattice in Figure [7J As we can see, these are just fluctuating about 
constant average values. Therefore, as with the JSD calculation, we conclude 
that the flip time distribution obtained during the PBC calibration stage is already 
very close to the steady-state flip time distribution, so our simulations became 
self-consistent very shortly after SBC is turned on. 
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Figure 7: The first to fourth moments of the T-to-| flip time distributions in the SBC0 simulation 
at T = 3.0 for a 100 x 100 square Ising lattice, versus the number of 10,000-Monte-Carlo-step 
sweeps. 

3.2. Comparison Against PBC 

To compare SBC0 and SBC1 against PBC, we simulate each boundary con- 
dition and each temperature ten times. For each simulation, the specific heat and 
magnetic susceptibility were each calculated from 10,000 independent data points. 
Since we have the analytical expression for the specific heat of a square Ising lat- 
tice |0], we divided the numerical specific heats by the analytical specific heat, and 
plot the ratio in Figure [8] As we can see, the SBC0 and SBC1 specific heats are 
slightly smaller than the analytic specific heat, while the PBC specific heat fluc- 
tuates about the analytical specific heat. Comparing the two SBCs, we find the 
SBC1 specific heat closer to the analytic specific heat. However, these differences 
are not strongly significant in statistical terms. 

In Figure |9] we show the magnetic susceptibilities x obtained for the three 
boundary conditions, divided by the magnetic susceptibility obtained using the 
high-temperature series expansion 16511 . Here we find that the SBC0 and SBC1 
magnetic susceptibilities are both lower than the high-temperature series expan- 
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Figure 8: Normalized specific heats measured from PBC, SBCO and SBC1 simulations of a 100 x 
100 square Ising lattice. We ran 10 simulations for each boundary condition and each temperature. 
For each simulation, the specific heat c is calculated from 10,000 statistically independent values 
of the average energy per spin e sampled every 2t. The error bar for c is then the standard deviation 
over the 10 equivalent simulations. 



sion, while the PBC magnetic susceptibility is slightly higher. As with the specific 
heat, the SBC1 magnetic susceptibility is closer to the series expansion magnetic 
susceptibility. None of these differences are statistically significant. Finally, we 
show in Figure HO] the spin- spin correlation G(r) obtained at three different tem- 
peratures for the three boundary conditions. For all separations, G(r) for the three 
boundary conditions are identical to each other to within numerical uncertainties. 

Through these measurements, we see that the SBC is just as good as the PBC 
for simulating the two-dimensional square Ising lattice in the high-temperature 
regime. In the next section, we will illustrate the advantages SBC have over PBC 
for Ising simulations, by considering Ising lattices with one special boundary. 

4. Asymmetric Lattices 

One very powerful application of SBC is that it can be used to simulate asym- 
metric lattices, for example, the Ising lattice in a spatially varying magnetic field. 
In these lattices, translational symmetry is broken, and PBC can no longer be im- 
posed on all four boundaries. However, so long as the local flip time statistics 
do not vary too quickly from point to point in the Ising lattice, we can still flip 
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Figure 9: Normalized magnetic susceptibilities measured from PBC, SBCO and SBC1 simulations 
of a 100 x 100 square Ising lattice. We ran 10 simulations for each boundary condition and 
each temperature. For each simulation, the magnetic susceptibility x is calculated from 10,000 
statistically independent values of the average magnetization per spin m sampled every 2t. The 
error bar for^ is then the standard deviation over the 10 equivalent simulations. 




Figure 10: Spin-spin correlations measured from PBC, SBCO and SBC1 simulations of a 100 x 100 
square Ising lattice at T — 3.0,4.0,5.0. For each boundary condition and each temperature, we 
repeat the simulation ten times to estimate the error bars. 
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pseudospins using the flip time statistics of systems spins on the boundary. 

In this paper, we investigated the Ising lattice with one boundary coupled to 
static t spins, as well as one with one open boundary. We call these the magnetized 
boundary (M) and open boundary (O) respectively. For both lattices, we expect 
the influence of the special boundary to be strong close to it, and weak further 
from it. Because spins in the same row are equidistant from the special boundary, 
we expect flip time statistics to be identical within a row, and slowly varying from 
row to row. This means that instead of two flip time distributions for SBCO or four 
flip time distributions for SBC1, we will need to work with a larger number of flip 
time distributions, at different distances from the special boundary. 




(a) first calibration stage (b) second calibration stage (c) implementation stage 



Figure 1 1 : To simulate an asymmetric Ising lattice with magnetized or open boundary conditions, 
we need to run the first calibration stage with PBC on all boundaries. Once we collect sufficient 
spin flip statistics, we switch over to the second calibration stage, where we impose the magne- 
tized/open boundary conditions on one boundary, SBC on the opposite boundary, and PBC on 
the remaining two boundaries. In this second calibration stage, we collect spin flip statistics row 
by row, and also distinguish between spin pairs parallel or perpendicular to the special bound- 
ary. In the implementation stage, we turn on SBCs that vary with distance away from the special 
boundary. 

4.1. Algorithm 

Two calibration stages are needed for asymmetric lattices, as shown in Fig- 
ure [TU In the first calibration stage, we collect the system-wide spin flip statistics 
within the lattice with PBC imposed on all four boundaries. After the system- wide 
spin flip statistics is obtained, we impose magnetized/open boundary conditions 
to one boundary and SBC to the opposite boundary. We continue to impose PBC 
on the other two boundaries. In this second calibration stage, we collect spin flip 
statistics row by row, to allow them to vary with distance from the special bound- 
ary. In this calibration stage, we also distinguish between spin pairs parallel or 
perpendicular to the special boundary, for higher-order SBCs. Because we now 
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need to construct many more flip time distributions, we use many queues contain- 
ing 100,000 flip times instead of one queue containing 1,000,000 flip times. 

Finally, in the implementation stage, we implement SBCj_(l) one row from the 
special boundary, SBC ± (2) two rows from the special boundary, and so on and so 
forth till SBCj_(L) on the row of system spins furthest from the special boundary. 
For the row of pseudospins along the opposite boundary, we impose SBC||(L). We 
allow the simulation to become self-consistent before performing measurements. 




Figure 12: Row magnetization per spin ;«(/) as a function of the distance r from the magnetized 
boundary, for a 100 x 100 Ising lattice with a magnetized boundary simulated using SBC0 and 
SBC1 at T = 3.0,4.0,5.0. 

4.2. Comparison Between Orders 

For the magnetized boundary and open boundary simulations, there are no 
analytical results to compare against. Therefore, we first compare the results for 
SBC1 against that of SBC0 for the Ising lattice with a magnetized boundary. In 
Figure [121 we show the row magnetization per spin, measured using SBC0 and 
SBC1. As we can see, the results of the two SBCs agree very well with each 
other. 

We also compare the transverse spin-spin correlation function measured us- 
ing SBC0 and SBC1 for an Ising lattice with a magnetized boundary, and also 
against the SBC0 spin-spin correlation function for a symmetric lattice with no 
special boundaries. This is shown in Figure [T3] As we can see, G ± (r) for SBC1 
agrees very well with that for SBC0. We also see that G ± (r) in the middle of the 
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Ising lattice with a magnetized boundary coincides with G(r) of the symmetric 
Ising lattice with no special boundaries. This is expected, since the middle of the 
asymmetric lattice is far away from the magnetized boundary. 
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Figure 13: The transverse spin-spin correlation G ± (r) in the top, middle and bottom rows for a 
100 x 100 Ising lattice with a magnetized boundary above the top row. In this figure, we also 
show Gj_(r) for the symmetric lattice with no special boundaries. For r > 6 our results are not 
reliable, because we used short queues containing only 100,000 flip times to implement the SBCs, 
as opposed to 1,000,000 flip times in the queue for the symmetric lattice. 

More interestingly, we find that Gj_(r) is suppressed in the row adjacent to the 
magnetized boundary, and also the row adjacent to the opposite stochastic bound- 
ary. Suppression of G ± (r) right next to the magnetized boundary is expected, 
because we need to subtract from (SjSj) the product (which is larger 

close to the magnetized boundary). The suppression of Gj_(r) in the row adjacent 
to the stochastic boundary gives us a sense of how well the SBCs are working. In 
principle, if we use an SBC with a high enough order, G ± (r) in this row would be 
equal to G ± (r) measured in the middle of the lattice, i.e. as if no boundaries were 
present. Indeed, the SBC1 G ± (r) for this row is closer to the bulk G ± (r) than the 
SBCO G_i_(r), giving us confidence that this convergence is indeed happening. 

For other physical quantities, the SBCO values also agree very well with their 
SBC1 counterparts. Therefore, from this point on, we use only the SBC1 results 
to compare the effects of different boundary conditions. 
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Figure 14: Row magnetic susceptibility per spin as a function of distance from the boundary, mea- 
sured using SBC1 from a 100 x 100 Ising lattice at with no special boundaries (blue), a magnetized 
boundary (red), and a open boundary (black) at T = 3.0, 4.0, 5.0. 

4.3. Comparison Between Boundary Conditions 

In Figured!] we compare the row magnetic susceptibilities per spin at differ- 
ent distances from the magnetized and open boundaries, measured using SBC1, 
and also compared against the row magnetic susceptibility per spin at different 
distances from a given boundary for the SBC1 symmetric Ising lattice. As we 
can see, the row magnetic susceptibilities per spin for all three boundary con- 
ditions reach the same bulk value by a distance r = 9 from the boundary, for 
T = 3.0, 4.0, 5.0. Close to the boundary, we find the row magnetic susceptibility 
per spin of the SBC1 symmetric lattice dip below the bulk value by about 10%. 
Since we expect an infinite-order SBC to completely eliminate the presence of a 
boundary, this gives a measure how 'imperfect' SBC1 is. 

For the Ising lattice with an open boundary, the row magnetic susceptibility 
per spin is suppressed to a level below the symmetric lattice. Based on the error 
bars found in our measurements, this suppression of the magnetic susceptibility 
by the open boundary is statistically significant. In particular, for T = 5.0, where 
SBC1 is effectively as good as SBC-oo, we find no visible suppression of the row 
magnetic susceptibility per spin at the boundary of the symmetric lattice. The row 
magnetic susceptibility per spin of the asymmetric lattice with an open boundary, 
on the other hand, is clearly suppressed at the boundary. The strongest suppression 
of the row magnetic susceptibility per spin occurs for the asymmetric lattice with 
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a magnetized boundary. This suppression is especially strong at T = 3.0, which is 
very close to the critical temperature T c = 2.269 of the square Ising lattice |45|]. 




Figure 15: Longitudinal spin-spin correlation G\\(r) between spins in the row adjacent to the 
boundary, and spins r + 1 rows away for a 100 x 100 Ising lattice at T — 3.0,4.0,5.0. In this 
figure, G]|(r) of the symmetric lattice with no special boundaries are shown in blue, G\\(j) of the 
asymmetric lattice with a magnetized boundary are shown in red, while G\\(r) of the asymmetric 
lattice with an open boundary are shown in black. 

Next, we show in Figure [15] the longitudinal spin-spin correlation Gy(r) be- 
tween spins in the row adjacent to the boundary, and spins at various distances 
from the boundary. Here we see that at higher temperatures (T = 4.0 and T = 5.0), 
Gy(r) coincides for all three boundary conditions. At T = 3.0, which is closest to 
Tc, G\\(f) for the asymmetric lattice with open boundary and G(r) for the symmet- 
ric lattice coincides with each other, while Gy(r) for the asymmetric lattice with 
magnetized boundary is suppressed. 

Finally, we show in Figure \\6\ the transverse spin-spin correlation function 
G ± (r) between spins at the same distance r from the boundary, for a lOOx 100 Ising 
lattice at T = 3.0. Here we see that G ± (r) measured in the middle of the lattice 
converge to G(r) of the symmetric lattice for both magnetized and open boundary 
conditions, i.e. the effects of the boundary conditions, whatever they are, are 
negligible this far away. For the row of spins adjacent to the stochastic boundary 
opposite the special boundary, G ± (r) is suppressed by the same extent for both 
boundary conditions. Since these spins are all adjacent to SBC1 pseudospins, we 
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judged that this common suppression is due to the 'imperfect' SBC1 algorithm. 
More importantly, for the row of spins adjacent to the special boundary, G ± (r) is 
suppressed beyond the SBC1 imperfection only for the asymmetric lattice with a 
magnetized boundary. 




Figure 16: Transverse spin-spin correlation Gj_(r) between spins the same distance r from the 
boundary of a 100 x 100 Ising lattice at T — 3.0. In this figure, G(r) of the symmetric lattice 
with no special boundaries is shown in green, Gj_(r) of the asymmetric lattice with a magnetized 
boundary are shown in red, while G±(r) of the asymmetric lattice with an open boundary are 
shown in blue. 

5. Conclusions 

To summarize, we introduced in this paper a hierarchy of stochastic bound- 
ary conditions (SBCs) for Monte Carlo simulations of the two dimensional Ising 
model. The main idea behind our SBCs is to couple spins at the boundaries of the 
system to independent pseudospins, whose dynamics and local correlations mim- 
ick spins in the bulk of the system. We then described the zeroth (SBC0) and first 
(SBC1) order algorithms in detail, in particular how we collect spin flip statistics 
in a PBC calibration stage, before turning on SBC in the implementation stage, 
and allowing the simulation to achieve self-consistency before measurements are 
performed. At simulation temperatures above T c , we find that self-consistency 
is reached very quickly. From measurements of the specific heat, magnetic sus- 
ceptibility, and spin-spin correlations, we found the two SBCs compare favorably 
against the PBC, and that SBC1 systematically improves upon SBC0. 
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We then demonstrate the advantage of using SBCs to simulate two asymmet- 
ric Ising lattices, one with a magnetized boundary, and another with an open 
boundary. To do this, we allow the SBC to vary with distance from the special 
boundary. We then checked from measurements of row magnetization per spin 
and transverse spin-spin correlations that the modified SBCO and modified SBC1 
agree with each other, and also with the bulk results expected far from the special 
boundary. To understand the effects of the magnetized and open boundaries, we 
measure the row magnetic susceptibility per spin, the longitudinal spin-spin corre- 
lations, and the transverse spin-spin correlations. We find that the magnetized and 
open boundaries both suppress the row magnetic susceptibility near them, with 
the magnetized boundary more strongly so. For the longitudinal and transverse 
spin-spin correlations, only the suppression near the magnetized boundary can be 
clearly seen in the SBC simulation results. 
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